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I. INTRODUCTION 

The research field of plasmonics^ has seen an enor- 
mous development over the last decades, both exper- 
imentally and theoretically. Examples of the aspects 
studied include enhanced light emission under vari- 
ous circumstances;^ enhanced spectroscopies such as 
surface-enhanced Raman scattering (SERS),— extraordi- 
nary transmission of light^, and biosensing applications.^ 

A corresponding development has also taken place on 
the theory side and a variety of methods are used to solve 
theoretical problems in plasmonics. These include exact 
methods that apply for certain geometries such as Mie 
thcorySr— for problems with spherical symmetry, but in 
most situations methods that make more extensive use of 
numerical calculations are needed. The finite-difference 
in time-domain method (FDTD}i2i"— is one such method 
that has grown in popularity in recent years, the discrete- 
dipole approximation (DDA) metho d^"^'^^ and Green's 
function (GF) methodi^"— are two other methods that 
are often used. Of these the DDA method has a some- 
what longer history and probably a bigger user base. The 
Green's function method on the other hand can be more 
flexible in certain situations. 

In this paper we will present a calculation of the 
Green's function for a layered material. In particular 
this makes it possible to study scattering off embedded 
inclusions such as nanoholes in metal films. ^^"^^ Paulus 
et al. presented a derivation of the Green's function 
(Green's tensor) for a layered system in Ref. [3. The 
present derivation follows the same basic ideas, but we 
derive rather elegant, explicit expressions for the Green's 
function that only involve a single transfer matrix recur- 
sion relation, and which makes it possible to explicitly 
demonstrate various symmetry properties of the Green's 
function such as reciprocity symmetry. We also study 
the analytic properties of the Green's function in Fourier 
space and show how this effects the long-range properties 
of the Green's function which for metallic films are domi- 
nated by plasmon polaritons for distances typically in the 
range of 100 nm to 10 /xm and for even larger distance the 



dominating contribution comes from a boundary wave. 

We will illustrate the Green's function method by cal- 
culating scattering cross sections for light off nanohole 
systems in thin metal films. The nanoholes of these sys- 
tems typically have diameters that range from 50 to 100 
nm in a thin Au film of thickness 20 nm. The opti- 
cal properties of nanohole system have attracted intense 
interest for over a decade now in the context of extraor- 
dinary transmission through an array of nanoholes dis- 
covered by Ebbesen and coworkers;^ and studied by var- 
ious theoretical methodsj^i^ for a couple of recent re- 
views on this subject see Refs. [s^ andlssl. But nanohole 
systems are also studied in connection with biosensing 
applications, since they at the same time can act as 
capturing centers for biomolecules and light scatterers 
whose properties are modulated by the presence of these 
molecules. The basic optical scattering properties of in- 
dividual nanoholes and chains of nanoholes in thin metal 
(Au films) have been studied by Kali and co-workers and 
the results show signs of strong hole-hole interactions 
Here we present theoretical results for the scattering cross 
section off multi-hole systems that are in good agreement 
with the experimental ones. The theoretical results com- 
bined with an analysis of the behavior of the Green's 
function shows that the hole-hole interaction affecting the 
light scattering is to a large extent mediated by surface 
(interface) plasmons. The special nature of the plasmons 
means that there is a strong relation between polariza- 
tion and propagation direction; hole-hole interactions are 
much stronger in the case when the electric field is polar- 
ized along the axis through the hole centers than when 
the polarization is perpendicular to the chain axis. 

The rest of the paper is organized in the following 
way. In section |TT] we give a brief overview of the Green's 
function method. Section IIIII details our calculation of 
the Green's function for a layered background through a 
transfer-matrix method and we also show how a number 
of physical quantities can be derived from the Green's 
function. Section IIVI focuses on the analytic properties 
of the Green's function in wave vector space and its con- 
sequences for the long range behavior in real space. Sec- 
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tion fVl gives a brief description of the numerical solution 
of the integral equation determining the electric field in 
the scatterers. In Sec.|Vl]we apply the method to a study 
of the optical properties of nanoholes in a thin metal film, 
and the paper is summarized in Sec. IVIII 



II. BASIC TREATMENT OF THE 
SCATTERING PROBLEM 

We consider first a situation where all of space is filled 
with a material with dielectric function eb, correspond- 
ing to a wave number 



(1) 



where fco and c are the wave number and speed of light 
in vacuum, respectively, and cj the angular frequency 
of electric and magnetic fields. The task at hand is to 
solve Maxwell's equations which assuming all fields have 
a e~*'^* time dependence, read 



along with 



V X £; = iujB, 



V X B = fioj - iujfioeoSBE, 



V • B = 0, and W ■ D 



(2) 



(3) 



(4) 



The solution for the electric field can in this case be 
written as the sum of a source term that depends on 
the current at the field point r and another term that 
through a Green's function takes into account the effects 
of the currents everywhere else^ 



where R = f — r' , and (g) denotes a dyadic product. 

In case the dielectric function is not constant in space, 
there is a modification of the second of the Maxwell's 
equations 

V X HoJ i^E ^ '-E, (9) 

c 

where the last term is new compared with the case of a 
homogeneous medium, and Srei can vary in space in an 
essentially arbitrary way. We rewrite this as 



V X B = ^oJtot 



(10) 



where the total current due to both external sources and 
scattering is 

Jtot = J + Jscatt = J 5 E. (11) 

Mo 



By letting jtot take the place of j in the solution we 



get 



E{r)= L-E(r)+ Gh[r,r) 



X iiiofioiif') + k^Ae{r')E{r')) d^r' 



(12) 



We will next show how this expression, in particular the 
Green's function, can be generalized to deal with a lay- 
ered system where the background dielectric function 
ssir) varies stepwise along one direction (z) in space, 
and Ae = erci — £b then describes further variations of 
the dielectric function due to scatterers. 



E{f) 



L ■ /(r) 
iueoSB 



G,.(r-,r').j(r')dV. (5) 



Vj-Vs 



III. GREEN'S FUNCTION FOR A LAYERED 
STRUCTURE 



Here L is a tensor which depends on the shape of the 
excluded volume Vs- In the most common cases where 
the excluded volume is cubic or spherical in shape L is 
diagonal and each of its three elements have the value 
1/3. The Green's function is given by 



VV 

Kg 



G^(f,r') = 



VV 

Kg 



47r|r — r'l 
(6) 

where is the Green's function to the scalar Helmholtz 
equation, and thus satisfies 



(V' + fc|)G^(r,P) = -5(3)(^_^)_ 

More explicitly we have 

ikBR — 1" 



(7) 



Gh{R)= 1^1 

3 — iikBR 



k%R^ 



1 



fr2 R2 



R 



AttR ' 



(8) 



A. Formulation in terms of 2D Fourier transform 

The Green's function for a homogeneous background 
given above in Eq. ([5]) can be written in terms of a Fourier 
integral as 

^'^^^ - J (2.)3fc| q^-kl ' ■ 

This expression becomes more useful in handling lay- 
ered structures if we integrate out the qz variable, which 
yieldsiS 

G.(i?) = -^'5(^ni?)+/(giG.(<zi|,z)e^?rn,. (m) 
<-> ^ ^ 

The notation Gh indicates that we are still dealing with 
the Green's function for a homogeneous background, but 
the formal generalization of Eq. (fH)) to G valid for a 
layered background is straightforward. The (5-function 
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term is the result of a subtract-add operation necessary 
to render the contour integral convergent. However, in 
the followin^we will not explicitly deal with the singular 
behavior of G when r — >■ f*, so we will therefore leave out 
this term in the rest of the calculation. The 2D Fourier 
transform (FT) of the Green's function in Eq. ([T4| is 



Gh{q\\,z,z') = 



2p 



1 - 



P 



(15) 



where p stands for the absolute value of the z component 
of the wave vector q (originating from the residue at the 
pole in the contour integration) 



P 



(16) 



The square root function in Eq. (|16|) should, to give phys- 
ical results in the form of outgoing, damped waves, be 
evaluated with the branch cut along the positive real axis 
of the argument. The superscript r on the wave vector 
q"^ indicates the direction of propagation of the waves 
(called primary propagation direction in the following) 
r = +1, or just r = +, when z > z' and r ~ —1 when 
z < z' . For the wave vectors we have 

= q\\ ±pz. = {{q\\/kB) cos0g, ((?||/fcs) sin(^,, zbp/fcs), 

(17) 

where q\\ = \q\\\- The corresponding unit vector, 

r = r/kB, (18) 

together with the unit polarization vectors for s polariza- 
tion 



\z X q 
and p polarization 



±1 



(-sin cos 0g,O), 



(19) 



P 



X (p^ 



±— cos (pg,±— sm (pq, - — 
Kb Kb Kb 



(20) 

form an orthonormal basis. The unit tensor therefore can 
be written 



1 = (g) g^ + p"" (g) + (g) 



(21) 



Using Eqs. (|18l) and I^T^ we can rewrite the Green's func- 
tion Fourier transform in Eq. (1151) as 



Gh{q\\,z,z') + f 

2p 



(22) 



As a first step towards generalizing Eq. (|22|) to a sit- 
uation with a layered background, we conclude that a 
particular element of the tensor can be written 



Gh,ap{q\\,z,z') =a- E = 

= a ■ \{p''A''P + srj^T,syTp{z^z') 



(23) 
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FIG. 1. Illustration of the plane waves generated in the dif- 
ferent layers of a four-layer system when the source is placed 
at z = z' . 



where £ is a vector field proportional to an electric field 
generated by the source. The wave amplitudes are found 
by projecting the source unit vector /3 onto the p and s 
unit vectors, thus 



2p 



p'' ■ P and A'^ 



2p 



13. 



(24) 



B. Generalization to a layered material 

When we turn to a layered material the source will 
still generate outgoing plane waves just like the expres- 
sion in Eq. (j23p indicates, however, now there will also 
be other waves reflected and transmitted at the different 
interfaces. 

This is illustrated in Fig. [1] for a system with four 
layers and the source placed in layer 3. In layer 3, there 
are waves going upwards and downwards on both sides of 
the source, in layer 2 there are also waves propagating in 
both directions, but in the two outermost layers there are 
only outgoing waves, propagating upwards in layer 1 and 
downwards in layer 4. The vector field £ we introduced 
in Eq. (1^51) takes the generalized form 



{pJAj^P{zoj)+sAj^\zo,i))e^^P'^^'^o,) 



+ {pJB;'^{zo,i)+sB;-'{zo,i)) e 



-irpi (z-zo,i) 



, (25) 



in layer I. Here r is the opposite direction of r, — r, and 
the unit vectors for p polarization as well as p depends 
on the layer number / since the background wave vector 
magnitude fcs varies from layer to layer. The wave am- 
plitudes AJ''^{zq i) for outgoing waves propagating away 
from z — z' and BJ''^{zq^i) for returning waves propa- 
gating towards z = z' in Eq. (j25p depend on the layer 
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number I, primary propagation direction r, and polar- 
ization a (p or s). The offset points Zqj are local origins 
for the plane wave exponentials, which can be moved 
around provided, of course, that the wave amplitudes 
are adjusted accordingly. The standard choice for zqj, in 
particular in a numerical implementation, is to use the 
bottom of all the layers above the source, and the top of 
all the layers below the source. In Fig. [T]this means that 
zo,! = di, zo,2 = d2, zo,3 = z', and 2:0,4 = d^. 

We need to determine the wave amplitudes Aj'"{zo^i) 
and Bj''^ {zq^i), something we will do in two steps. First 
we view the stack of layers as made up of two indepen- 
dent parts, one above the source plane z = z', and one 
below. We introduce relative wave amplitudes a[''^(zo_;) 
and bJ''^{zQ^i), corresponding to the actual amplitudes 
Aj'"{zoj) and Bj'"^ (zq^i). The actual amplitudes above 
the source are found from the relative ones as 



a{,{z'^ 



a{,[z') 

(26) 

where V denotes the source layer and z" is any z coordi- 
nate. In the same way, below the source 



ATiz") 



ai{z") 



&r(^") 



a;, [z') 



-ArAz'), B7{z") = ' ^ ' AtAz') 



(27) 

The relative amplitudes can be determined by a transfer- 
matrix calculation using the fact that there are only out- 
going waves in the outermost layers. But Eqs. (p6)) and 
([27]) show that the actual wave amplitudes and Aj, in 
the source layer play the role of driving forces for all the 
waves above and below the source, respectively, and still 
have to be calculated independently. This is done in the 
second step of our calculation, by a detailed investigation 
of the situation in the source layer. 



C. Transfer matrix calculation 

The calculation of the relative amplitudes uses the fact 
that there are no returning waves in the outermost layers, 
thus 



B+ = 1)1 = 0, and B^ = b^= 0. 



(28) 



We set the relative amplitudes of the outgoing waves in 
the outermost layers, i.e. al and a^, to 1, 



dl{zj^) = 1, and aj^{z-) = 1. 



(29) 



The coordinates z+ and z_ lie above {z^) and below 
(z_) all interfaces as well as the source and field points, 
respectively. 

The remaining relative amplitudes can then be deter- 
mined by recursion, applying the Fresnel formula at each 
interface and adjusting the amplitudes by exponentials to 
account for the propagation through intermediate layers. 



For a general z coordinate z" we get 



iz") 



— T/f/+ 



W+iz",z+) 



S+{z",di-,)T+^^ 



above the source, and 



"-1.1-1 



■T+fS+{d,,z+) 



(30) 



a 



bi-'^" 



iz") 

Si{z",di)T{'i:^^ 



W-{z\z.) 



■T, 



N~l,N'-'N 



Sff{dN-l, Z-) 



(31) 



below the source. Here is a "total" transfer matrix 
built up by factors 5, related to the wave propagation 
in the different layers, and T, describing reflection and 
transmission at a particular interface. 

The propagation in one single layer just yields expo- 
nential factors multiplying the wave amplitudes. We have 



a 



bi'^iz) 



SJ{z,z") 



al%z") 
bl'^z") 



where a denotes a polarization (s or p) and where 



Sf{z,z") = 









= TiPi(^-z") 



(32) 



(33) 



The cross-interface transfer matrices T relate the wave 
coefficients on opposite sides of an interface z = d, sepa- 
rating layers I and n (where Z = n ± 1), to each other 



(34) 



They can be evaluated by using the Fresnel formulae for 
s and p polarized waves. We express the result in terms 
of the reflection amplitudes for s and p polarized waves, 
respectively, incident from the material in layer I onto 
material n in case this is the only interface. 



aj'^id) 




' a]f{d)' 


. b^id) 




bl^-{d) 



fin — 



Pi - Pn 
Pi +Pn 



and 



SnPl ~ £lPn 
£nPl + £lPn ' 



(35) 



These quantities depend on the dielectric functions £/ 
and En and wave vector z components pi and Pn of the 
two layers. The T matrices are found after some alge- 
bra which also involves the amplitude of the transmitted 
wave. We get 



rpTS 

^In — 



1 



1 fin 

fL 1 



for s polarization, and 



rpTp klPn 

^in ^ 



1 



knPl 1 - fL 



1 fP 

^ Jin 

fP 1 

Jin ^ 



(36) 



(37) 



for p polarization. 
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D. Calculation of the primary wave 



we can write (using as a general polarization vector, s 
or p) each of these amplitudes as 



Using the scheme outlined in Sec. IIII C] we can calcu- 



A[f{z')^—al-P + x'^%z')Al-%z'), r = ±l. (40) 



late all the relative wave amplitudes we need. Now it 

remains to find the primary wave amplitudes A^'" in the This system of two equations has the solution 
source layer. We have already seen that in the case of a 

homogeneous material we have ^ ^ (^^rj, + aj, x^^" {z')) ■ fi i 

^ - l-x+^^{z')x-^^{z') 2^ ^^^^ 

in which the first term in the numerator is the direct 
wave from the source, the second term is the wave re- 
fiected once off the opposite interface, and the denomina- 
tor accounts for repeated reflections off the surrounding 
interfaces. We can now calculate the Green's function by 
using this solution in Eqs. ([25]), ([Ml), and (l27l) . 



A;;P{z')^—P^-$ and Aji' {z') ^ — ■ p. {38) 



In the present case we have to add the wave reflected 
off the "opposite" interface of the source layer to each 
of these expressions. Thus, the amplitude of the wave 
propagating upwards has a direct contribution from the 
source, and one contribution from the interface below the 
source, and vice versa for the wave propagating down- 
wards. By introducing the response functions, i.e. the 
ratios between reflected and incident wave amplitudes, 



E. Results for the Fourier-space GF 

With a source pointing in the /3 direction we can 
now write down the result for the matrix element 
Gap {q\\,z,z') = a- £ in Fourier space. To keep the whole 
thing manageable we divide the Green's function into one 
p part and one s part 



G(q|| , z, z') = G-P{q\\ , z, z') + G'{q\\ , z, z'), 



(42) 



X^'"(^") 



&[--(z") ^ i3['-(z") 
ar(z") Ar(z") 



(39) 



with 



GP( 



,z,z') = 



and 



G'{q\\,z,z') 



' a[;P(z') aJ;P{z') 



al^z') + a[;^(z') 



[pl + X^'nz')pl) « 



l-X+'P{z')x-'''{z') 2pi 
[l + X~''{z')] 



TT- (43) 



l-X+''{z')x-^'iz') 2pi 



TT-- (44) 



These expressions are a good starting point for a numerical implementation. 

For theoretical purposes it is, however, quite useful to express the Green's functions in terms of the W transfer 
matrices. By multiplying the numerators and denominators in Eqs. (|43p and (|44|) by aj/"^ {z')aj,''^ {z'), and using that. 



in view of Eqs. ([50)1 and pip . aj,f{z") and bj,f{z") are the flrst column elements of W'^'"'{z", Zt) we arrive at 



G?'(gi|,z,z') 



W^fiz,Zr) [p[ +p[x"-n^)] ^ [PI' +PlX^^nz')] W^fiz',Z-) 

D^z') 



and 



where 



,z,z') 



W^i^z, Zr) [1 + X^'^jz)] s^§[l + x^^'jz')] Wl{\z', Z-) 
Dfiiz') 



Dl,{z') = ~2ipi> [wt{%z',z+)W^^^{z',z^) - W+{%z',z+)W^f{z',z^)\ 

I 



(45) 



(46) 



(47) 



is the 1 , 1 element of the matrix 

1 1 



D"{z')^ [W+^^ {z' ,z+)Y{-2ipi,a,)W-^" {z' ,z^). (48) 



Here the superscript t denotes matrix transposition, and 
Gz is the Pauli matrix 



1 

-1 



(49) 
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To summarize, Eqs. p5)) -(l47 |) show how the Green's func- 
tion can be uniquely expressed in terms of the transfer 
matrices describing propagation from the outer layers of 
the system to the source and field points z' and z, re- 
spectively. 



F. Calculation of the real-space GF 

The layered system we are considering has cylindrical 
symmetry and it is therefore quite natural to view both 
the Green's function in real space, as well as its Fourier 
transform which has been at the center of our attention 
so far, as functions of cylindrical coordinates 

G{r, EE G(p, z, z') and G[q\\ , z, z') = G{q\\ ,<j)q,z, z'), 

(50) 

respectively. 

Thanks to the cylindrical symmetry the Fourier trans- 
form of the Green's function for an arbitrary (j)q can be 
related to the one a.t (pq = through 



G(g||,0„z,z') = U{cl>q)G{qn,cf>q = 0,z,z')[U{cbq)Y, 



where 



Ui(t>q) 



COS 4>q — sin 4>q 



sm( 




COS (j)q 







1 



(51) 



(52) 



and the superscript t denotes transposition. For (j)q = Q 

G^(g|| , 0, z, z') has four non-zero components xx, xz, zx, 

and zz, while for G'*(g|| , 0, z, z') only the yy component 
is non-zero. 

Likewise, for cj) — the Green's function in real space 
has 5 non-zero components xx, yy, xz, zx, and zz and 
for a general we have 



G{p, 0, z, z') = U{(f,)G{p, - 0, z, z')[um 



(53) 



Therefore to calculate G{p, 0, z, z ) in practice, we first 
calculate G for (f> — using the generalization of Eq. 
to the case of a layered material and then use Eq. ([55)) 
to get the final result. The angular part of the Fourier 
integral can be carried out analytically by making use of 
the integral representation of the Bessel functions 



i"Jn(^) 



1 



e*^'=°"'^cosn. 



(54) 



We get 

G(p,0,z,z') = 



(2^) 



■U{<Pq)G{q\\,0,z,z')[U{<t>q)Ye'^~rrii 



Gxxip,0,z,z') 


Gza:{p,0,Z,z') 



Gyy{p,0,Z,z') 





Gzz{p,0,z,z') 



(55) 



and using e^^r^ii = e^QnPcos^, ^ different 

components are explicitly given by 



Gxxip,0,z,z') 



i\\p 



xGxx{q\\,Q,z,z') + '^^-^^^Gyy{q\\,Q, z, z') 



q\\p 



dq\\ 
27r 



(56) 



Gxz{p,0,z,z')^ I iJi{q\\p)Gxz{q\\,0,z,z') —, (57) 



and 



GUp,Q,z,z')= I Jo(q||p)G,,(g||,0,z,z')^. (58) 



27r 

The expression for Gyy{p, 0, z, z') is obtained by the index 
replacements xx yy and yy — > xx in Eq. (j56p . and 
Gzx{p,Q, z, z') is obtained by replacing the index xz in 
Eq. ^ by zx. 

The integrations in Eqs. ([56l) . ([57]) and (|58p nominally 
run along the real g|| axis, however, as discussed in Ref. 
[Tsl the numerical evaluation can be speeded up by de- 
forming the integration contour into the complex plane. 

As a first step we also divide the Green's function into 

<-)• 

two parts, the homogeneous part Gh that we have al- 
ready discussed and given explicit expressions for, and 
an indirect part Gmd^ The homogeneous part refers to 
the situation without interfaces; the indirect part con- 
tains all contributions to the Greens function that at any 
point involves the reflection or transmission of a wave at 
any of the interfaces of the layered system. Both in real 
space and Fourier space 



G — Gh + Gi' 



(59) 



Thus, in practice we evaluate Gh from the explicit ex- 
pression in Eq. ([8]),'^- while the indirect part is calcu- 
lated from the Fourier integrals above. The functions 
in the respective integrals are analytic in the lower half 
plane (LHP) but has two branch cuts in the upper half 
plane (UHP) along the hyperbolas for which Lti[pi] — 
and Im[pAr] = 0, respectively, where pi and pn are given 
by Eq. (fT6|) using the material properties of the top and 
bottom layers. In addition the integrands may have one 
or several poles in the UHP. We deform the integration 
contour so that it starts from g|| = 0, first runs along 
the negative imaginary axis, then goes parallel to the 
real axis until it reaches a point beyond the singularities 
in the UHP where it goes back to the real axis. From 
there the integration either proceeds along the real axis 
to values of large enough that further contributions 
are negligible, or in the case that the lateral distance p 
between the source and field points is large a faster con- 
vergence is achieved by rewriting the Bessel functions in 
terms of Hankel functions asi^ 



1 r 



J„(<Zp)-- H^^\qp) + H(^\qp) 



(60) 
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and then carrying out the integration of the Hn part 
along a vertical path in the UHP and the Hn part along 
a vertical path in the LHP. 

For the purpose of (semi-) analytical calculations of the 
Green functions it is often an advantage to deform the 
branch cuts, more about that later. 



G. Reciprocity symmetry of the Green's function 

Reciprocity, which can be stated as "interchanging the 
source and the field probe does not change the result," is 



a central property of linear, time-reversal-invariant elec- 
trodynamics. In our case this requires that the Green's 
function fulfills the relation 

G/3a(r2,ri) = GQ/3(fi,r2). (61) 

To show that this is in fact true one can go back to Eqs. 
(|45l) and (j46|) . and first look at the matrix appearing 
in the denominators, D'^{z'). While D'^{z') nominally 
appears to be a function of the source coordinate z', it 
is in fact an invariant, independent of z' . To see this we 
first note that D'^ is independent of the position of z' 
within layer V . Writing D'^ as 



D'^iz') = [W^+'"(dp_i - 0, z+)] * S+{z', dv^i){-2ipvcj,)Si;{z\ dv)W-'''{dv + 0, z_), 
explicitly exposes the propagation in layer /' and it is easy to show that 

S'l^{z' ,dv-i){-2ipv(Jz)S^{z' ,dv) = [dv , dv -i)(~2ipv a z) = (-2ip/'cr^)S',7(d;'-i, d^). 
It remains to see what happens to D'^{z') when z' is moved across an interface. We then have 

D''idi,+0) = [W+^^{di, +0,z+)Y i-2ipi,a,)Ti^j,^^W-^^{di, -0,z^), 

just above an interface, and 

D''{di.-0) = [W+'''{di, +0,z+)]*Ti^+^^i,{-2ipi,+ia,)W-'^di, ^0,z_), 

just below. But for both p and s polarization an explicit calculation shows that 

I 



(62) 
(63) 
(64) 
(65) 



(-2ipi'Crz)TlTF+i = Ti,^^i,{-2ipi>+ia:,). 



(66) 



Thus, Eqs. (|63| and (|66|) show that the matrices and 
D" are invariant to all changes of z', both within a layer 
and from one layer to another. 

To prove Eq. (|5T|) we reverse the propagation direc- 
tion in Eqs. psj) and pS)) . which means that g|| — 7> — 
T — T, z and z' and the layer indices I and V are inter- 
changed, and s(g]|) — s(— (7||) = — s((7||), and p^{q\\) — 
p^(— (7||) — p*((7j|), and find that 



H. Surface response from the Green's function 

As we have seen in this section, calculating G for a lay- 
ered system is in general fairly involved, however, once it 
has been calculated a lot of information can also be ex- 
tracted from the Green's function and its Fourier trans- 
form. As a first example we determine the reflection fac- 
tors for a plane wave impinging on one of the outer in- 
terfaces at z = c?i or z = d^-i- 

For definiteness we concentrate on the reflection off the 
top interface at z = c?i, and use 



G/3q(^<?|| I ^2, Zi) — G^^((7|| , Zi, Z2) 



(67) 



Ga{i{q\\,z,z') = d • f , 



(69) 



for both p and s polarization. As a consequence, inserting 
Eq. (H^ into Eq. (ITil) and substituting the integration 
variable g|| — )■ —q\\ we recover Eq. 



(2^ 



■G/3c«(g||,Z2,zi)e*«i|-^'-^ii 



■('=■211-^11) 



(cf. Eq. ([23])) with a vector field £ of the form shown in 
Eq. ([^5]) to find a relation between the reflection coeffi- 
cients and the Fourier transform of the Green's function. 
We assume that the source is placed above the interface 
so that the vector field can be written in terms of the 
actual amplitudes we introduced earlier 



E = 



G^„(-q||,z2,zi)e''r(^^i||-^"^ii) = G„^(fi,f2), 

(68) 



Gf (gji , rfi + 0, z') + G" (9|| , rfi + 0, z') • /3 



(27r)2 

where we used Eq. ([57]) in the last step 



rA-^^P{d^) + p+B-^^\d^) + s-A-^\d,) + s+B-^\di). 

(70) 

Moreover, if the Green's function is divided into a ho- 
mogeneous part and an indirect part as in Eq. (j59p . in 
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this case the terms with A coefhcients contribute to Gh 
whereas the B terms contribute to Gmd- We can there- 
fore conclude that the reflection coefficient for polariza- 
tion a can be written 



Ka = — — = • [' ^) 

A, - (di) a- ■Gh{q\\,di+Q,z')-a- 



I. The Green's function and far-field calculations 

In a lot of situations one wants to calculate the scat- 
tered electric field very far from the layered system. 
Given a source distribution jtot(^) the field, retaining 
only non-zero terms in Eq. (jI2p . is 



E{r) = I G{f,r)iojfiojtot{r)d^r' 




10 10' 
Lateral distance p (nm) 



(72) 



FIG. 2. (color online) Behavior of the absolute square of the 
zz element of the Green's function along the surface of (a) a 
gold sample, (b) a 20 nm thick gold film. The photon energy 
is 1.8 eV. 



The Green's function here can be written 



G{r, 



■G((?i|,2;+,2: ) 



(27r)2 

xe*«i|-(*'i|-i)eV''Bi"«ii^'="^+^ (73) 



where z+ lies above all layer interfaces dj, as well as the 
source z' . To get somewhat simpler expressions we first 
assume that we can set z+ = 0. The integral can be 
evaluated by the method of stationary phase which yields 



G(f,r') = 



47rr 



(74) 



where 



^far{d, (p, z') = —2ikBi cos9G{kBi sin 6*, ip, 0, z'), (75) 

r = |r| and ksi — yfe\Lo j c. In case z+ > one must use 
the generalized expression 

gfar(0, V', z') ^ -2ikBi COS Oe-'P' '+G{kBi sin 9, if, z+,z'), 

(76) 

where the exponential function compensates for the prop- 
agation of the outgoing wave included in G. 

This expression for the far field is also useful in order 
to evaluate the field generated in the layered system by 
an incident plane wave. Assume that a plane transverse 



wave 



Ak-r 



impinges on the top surface of the stack of layers. This 
plane wave can be generated by a point source very far 
away at the point (r, 0, ip) in spherical coordinates in the 
direction where the wave comes from, i.e. 

k = (— fcsin6'cos(/3, — A:sin6'sin(/7, —kcosO). (77) 



Comparison with Eq. ([5]) shows that this requires a point 
source 



(78) 



at the point r. The full field at the point rg = (a^o, yo, zq) 
in the layered structure can now be calculated by insert- 
ing the source of Eq. ([78| in Eq. ([12]) (generalized to a 
layered background) and then applying the reciprocity 
relation Eq. (HI]), and Eq. ([71]) • This yields 



(79) 



IV. LONG-RANGE PROPERTIES OF THE 
GREEN'S FUNCTION 

The asymptotic behavior of the Green's function in the 
case that both the source and field points lie close to a 
metal surface is a problem that has attracted a lot of 
interest in the last few years j2§r^^ 

Our results for the amplitude squared of the zz element 
of the Green's function are shown in Fig. [5] (which uses 
a logarithmic scale on both axis). This figure illustrates 
that the field generated by a point source near a metal 
film and propagating outwards along the surface basically 
displays three different regimes: At short range from the 
source the dipole field originating from the source (and 
its image) dominates completely giving rise to a large 
field that however drops off as After that follows a 

range of distances where plasmon propagation along the 
surface gives the dominant contribution to the Green's 
function. The plasmons are cylindrical waves confined 
to the metal surface so their amplitudes decay as 1/ ^/p 
in the absence of power losses. A Au film on a glass 
substrate supports two different plasmons, and in this 
case we see interference between them. Eventually, once 
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FIG. 3. (color online) Deformation of the integration path in 
the complex plane in order to capture the asymptotic behavior 
of the Green's function along the surface of the metal film. 



we get to lateral distances between the source and field 
points corresponding to the plasmon propagation length 
the Green's function drops exponentially due to losses in 
the metal film and we enter the domain where the main 
contribution comes from a boundary wave, also known 
as a Norton wave^^ propagating along the interface. 

Naively one may expect IG^^p to decay as in 
the boundary wave regime, but in fact one finds a faster 
decay, ^ l/p'*- Refs. l36r,39. present a number of deriva- 
tions of this behavior. The basic physical reason behind 
the rapid decay is a destructive interference between the 
direct wave emerging from the source, and the wave re- 
flected off the surface. As is easily seen from the expres- 
sions in Eq. ([55]). exactly at grazing incidence (where 
pi = 0) both of the Fresnel reflection coefficients /" and 
equal -1, which means that to lowest order the sum of 
incident and reflected wave vanishes. Away from grazing 
incidence the incident and reflected wave do not cancel 
exactly and what remains (with |Gzzp ~ is the 

result of the interference between these contributions. 

In order to understand and calculate the asymp- 
totic behavior of the Green's function it J^s useful to 
study the behavior of the Fourier transform G{q\\ , 0, z, z') 
in the complex plane. Figure [3] shows the general 
structure for the case of a three-layer structure vac- 
uum/metal/dielectric. The FT has two branch points at 
ksi and fcsa, the vacuum and dielectric wave numbers. 
The physical branch cuts, as discussed in connection with 
Eq. (fTB)) . in this case follows the real axis back to the ori- 
gin and then runs out along the positive imaginary axis. 
For the purpose of performing the Fourier integral for 



large separations between the source and field points it 
is, however, better to deform the branch cuts to run par- 
allel with the imaginary as shown in Fig. [31 

Now the Bessel function in the integral Eq. (1551) can 
be split into two Hankel functions as indicated in Eq. 
([BUI) . The integral containing Hn [qp) is calculated by 
deforming the contour to run far into the LHP, which for 
large p yields negligible contributions. The integral con- 
taining H^i^ (qp) is calculated by deforming the contour 
to run far into the UHP where this Hankel function is ex- 
ponentially small so that the contributions from most of 
the contour are negligible. However, unlike in the LHP, 
the contour in the upper half plane must return to the 
real axis (or its vicinity) at every obstacle in the form of 
a branch cut or pole, and these parts of the integration 
path yield the dominating contributions to the Green's 
function for large The singularities closest to the real 
axis gives the contributions with the farthest range in p. 

For the case illustrated in Figs. [H and [31 the contri- 
butions from the integration along the branch cuts give 
the long-range, boundary wave contribution that persists 
over the entire range of distances in Fig. [H A close look 
at the result for a film shows that the long-range tail ex- 
hibits some small oscillations between the contributions 
from the two different boundary waves, the one at the 
vacuum-gold interface and the one at the gold-glass in- 
terface. 

As already said the boundary waves become the dom- 
inant contribution for values of p beyond the plasmon 
propagation length. In Fig. [2{b) this happens around 
p ^ 10 i-im. The behavior for p ~ 1-10 pm is dictated 
by two different plasmons corresponding to the two poles 
in the UHP. The pole to the right of both branch cuts 
corresponds to a charge-symmetric, bound plasmon with 
a wavelength of 364 nm, shorter than the wavelengths of 
huj = 1.8 eV light in both glass and vacuum."*^ This mode 
is thus bound to the film, evanescent in both vacuum and 
glass, and the pole lies on the physical sheet of the Rie- 
mann surface. The other plasmon with a wavelength 666 
nni is evanescent in vacuum, but propagating in glass. 
Therefore this charge-asymmetric mode is termed leaky 
since energy is lost to the glass side4^ From a formal 
point of view this mode does not correspond to a bound 
state and the corresponding pole lies on a higher sheet of 
the Riemann surface, but is brought out in the open by 
the contour deformation. 



V. SOLUTION OF THE SCATTERING 
PROBLEM 

The ultimate goal with the use of the Green's function 
method is of course to solve electrodynamics problems in 
the form of Eq. (|12p generalized to the situation with a 
layered background. 



L ■ Jjr) _ Asjr) 



L-E{r) 
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+ / G(r, r') (iujfiojir') + klAe{r')E{r')) d^r' . 

(80) 

We focus on the case where the current sources j are well 
separated from the scatterers so that the contribution 
from the current term in the integral in Eq. (jl2p can be 
written as an incident field, 

%{f)^l G{ry)iuj^i„j{r^)d\'. (81) 



If the sources are very far away we have a situation where 
a plane wave is incident on the layered system, and gives 
rise to reflected and transmitted waves in the system, all 
of this is dictated by the presence of the Green's function 
in Eq. ([5T|) . With an incident field written as EincG^^'^ 
we obtain a driving field Eo^r) in the layered system in 
accordance with Eq. (|79|) . 

Then Eq. can be written 



E{r)=Eo{r)-^^L-E[r)+ ( G{f,r^)BAe{r')E{r')d^r'. (82) 
esW Jv,-Vs 

By moving all terms involving the full field E{f) to the left-hand side, thus leaving only the driving field Eo^r) on 
the right hand side, and then discretizing the electric field on a mesh of equally sized cubic elements that covers all 
scatterers we arrive at 

-^qrs H ■ Eqrs ^o^^qrs^'^^qrs / ^ k^^Sq' r' s'^KlGq—q' _r — r' ,s,s' g' — E/Q^q^s. (§3) 

^B.qrs , , , 

11 q'r's' 



Here r, and 5 are discrete coordinates for the mesh 
elements in the x, y, and z directions, respectively. With 
a mesh side aM and an equivalent radius Rm the volume 
of a mesh element is 

VM^al,^^-^. (84) 

The term containing M describes self interaction within 
a mesh element. We use an approximation for M corre- 
sponding to a spherical mesh element of radius -Rm, 

M = [(1 - ikoRM) exp{ikoRM) - 1] t. (85) 

At the same time the self-interaction term is excluded 
from the sum (as indicated by the prime) since including 
this would involve singular contributions to the Green's 
function. 

Equation (1831) corresponds to a system of linear equa- 
tions; the left hand side can be seen as a 3Nm x 3Nm ma- 
trix multiplying a vector with 3Nm elements, Nm being 
the total number of mesh elements. We solve the system 
of equations iteratively using the stabilized biconjugate 
gradient method, BiCGstab(2).^'^ The iterative solution 
involves a large number of matrix multiplications. The 
contribution from the term in which the Green's func- 
tion multiplies the electric field is, as can be seen in Eq. 

the result of a convolution sum in the x and y di- 
rections. This means that the matrix multiplication can 
be speeded up by using a Fast Fourier transform (FFT) 
in these two directions'^ The same technique is used in 
the DDA methodJ' We calculate the Fourier transforms 



of the Green's function and the electrical field and mul- 
tiply the transforms by each other locally on the mesh 
in Fourier space and then transform the product back to 
the real space mesh. The use of the FFT is crucial in 
reducing computation times, since most of the compu- 
tational effort required in determining the electric field 
in the scattering volume goes into solving the equation 
system, thus essentially the repeated matrix multiplica- 
tions. The calculation of the Green's function, on the 
other hand, just needs to be done once per photon fre- 
quency, combination of z and z', and in-plane distance 
P- 

Once we have a converged solution to the system of 
equations the electric field inside the scatterers is known. 
At this point Eq. (j82p provides an explicit expression for 
the electric field everywhere else in space that can be 
evaluated by discretizing the integral as in Eq. (I55|) . 

The results presented in the next section focuses on the 
scattering cross section and thus depend on the far field 
which can be found from a discretizcd version of Eq. (j72p 
using Eqs. (fM)) and ([75]). The scattering cross section is 
given by 

do_ ^ r'^Sr 

d^ Sin 

where Sr is the radial component of the Foynting vector 
at a large distance r from the scatterers and Sm is the 
Foynting vector of the incident field. Given the (trans- 
verse) far field E{r), 

Sr = ^ceoV^\E{r)\\ (87) 
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where Eb is the dielectric function of the material the 
radiation is scattered into. 



VI. SCATTERING OFF NANOHOLES IN A 
THIN METAL FILM 

We now turn to calculating scattering spectra off 
nanoholes in a thin Au film. Such systems have been 
studied experimentally by Rindzevicius et al^ and 
Alaverdyan et ai— 



(a) 



E 




y (b) d 





1000 



FIG. 4. Illustration of illumination of two nanoholes with an 
electric field parallel to the dimer axis [in (a)], and perpendic- 
ular to the dimer axis [in (b)]. We also indicate how charges 
in the metal film surrounding the holes will be distributed 
in the case that the frequency of the incident light lies well 
below the single-hole resonance. One should note that the be- 
havior of nanoholes in terms of induced charges is essentially 
opposite to that of metallic nanoparticles. 

To study the problem theoretically we let a number of 
circular cylindrical holes in a Au film on top of a glass 
substrate act as scatterers. The Au film here has the 
same thickness, 20 nm, as in the experimental studies. 
The system is driven by a plane wave that impinges on 
the film (and the holes) at normal incidence, polarized 
either parallel to the symmetry axis of the hole chain 
or perpendicular to that symmetry axis, as illustrated in 
Fig. m We study primarily the forward scattering cross 
section as the edge-to-edge distance d, between the holes 
is varied. 

As a prelude we look at the Green's function with both 
the source and field points placed inside the metal film, 
which is a central quantity determining the interaction 
between different nanoholes in the film. Figure [5] dis- 
plays the behavior of the diagonal elements of G for the 
vacuum/ Au film/glass substrate system at a representa- 
tive photon energy of 1.8 eV (A ~ 690 nm) as a function 
of the lateral separation x. The xx element is by far the 
strongest over most of the range of distances x, between 
the source and field points. A source pointing in the x 
direction can excite plasmons propagating in the x direc- 
tion which explains why we have long-range interactions 
in this case. These plasmons are of the bound, charge- 
symmetric type discussed in Sec. IIV[ and illustrated in 
Fig. [21 The bound plasmon wavelength for hu — 1.8 eV 
is ~ 364 nm. The yy element is of comparable strength 
as Gxx for distances up to ~ 200 nm, i.e. in the near- 
field zone. However, for larger x the yy element is much 
smaller because a dipole pointing in the y direction can- 
not excite plasmons propagating in the x direction. Fi- 



FIG. 5. (color online) The diagonal components of the 
Green's function as a function of the lateral distance x along 
a 20 nm thick Au film on a glass substrate. Both the source 
and field points are placed in the middle of the film to best 
describe hole-hole interaction. The marks near the upper bor- 
der show where the 1st, 2nd, and 3rd neighbor hole is placed 
in a chain of holes with a radius of 40 nm and edge-to-edge 
distance d = 160 nm. (We stress, though, that the Green's 
function here has been calculated in the absence of any holes.) 



nally looking at Gzz we see that this component is much 
smaller than Gxx for sll x values. This is due to the 
boundary conditions for the electric field at a metal in- 
terface, which strongly suppress the normal component 
inside the metal. As a consequence of reciprocity this 
also means that a source inside the metal film oriented 
perpendicular to the interfaces is not very effective in 
generating electric fields elsewhere. 

The behavior of different elements of the Green's func- 
tion leads to differences in the hole-hole interaction de- 
pending on the polarization direction of the incident light 
(illustrated in Fig. 2]); interaction effects are much more 
important in the case of parallel polarization. The con- 
sequences are clearly seen in Figs. |6] and [T] Figure |6l to 
begin with, shows calculated scattering cross sections for 
two nanoholes of diameter 80 nm that are illuminated by 
light polarized parallel to the dimer axis. Each curve cor- 
responds to a different edge-to-edge separation between 
the holes. To make a comparison that brings out the ef- 
fect of hole-hole interactions the result for a single hole 
is also shown. This result is multiplied by 4 to adjust 
to the difference in scattering volume between the one- 
and two- hole cases. For a small separation between the 
holes the scattering cross section is suppressed and blue- 
shifted compared with the one-hole case. This is a re- 
sult of Gxx being negative for x smaller than « 200 nm 
(an edge-to-edge separation of 40 nm corresponds to a 
center-to-center distance of 120 nm). The shift can also 
be understood in view of Fig. HJa): the figure shows that 
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FIG. 6. (color online) Calculated forward scattering spectra 
for two nanoholes of diameter 80 nm in a 20 nm thick An film 
placed on a glass substrate illuminated at normal incidence by 
light polarized parallel to the dimer axis. The different curves 
show results for a series of edge-to-edge distances d between 
the holes as indicated in the key. The curve marked "1 hole" 
shows the corresponding result, adjusted for the scattering 
volume, for the case of a single hole. 
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FIG. 7. (color online) Forward scattering spectra for two 
nanoholes as in Fig. [S] however, here the incident light is 
polarized perpendicular to the dimer axis. 



for frequencies below resonance the field caused by the 
induced charges at one hole will counteract the exter- 
nal field at the other hole, but this situation is reversed 
for frequencies above the single-hole resonance hence the 
blue-shift. With an increasing distance between the holes 
the scattering cross section increases and its maximum 
red-shifts, and a maximum in the cross section occurs 



for d = 160 nm, corresponding to a distance of 240 nm 
between the hole centers. As can be seen in Fig. [5] 
this is close to the distance where Re[Ga;a;] has a max- 
imum. In this situation there is a constructive interfer- 
ence at one hole between the incident field and the field 
scattered off the other hole. For larger d the scattering 
cross section continues to red-shift, while the peak value 
falls off. For the largest separation d = 360 nm, we in 
fact see a new peak building up at the blue end of the 
spectrum (near 650 nm). For even larger d this peak 
grows and red-shifts reaching a second maximum around 
d = 560 nm corresponding to a center-to-center separa- 
tion right near the second maximum of Re[G2:a:] in Fig. [5] 
at X ~ 650 nm. In Ref. [H it is argued that the scattering 
from a chain of holes should show maxima whenever an 
odd number of half surface plasmon wavelengths can be 
fit in between two holes. We note that in the present 
case with Api = 364 nm, this predicts scattering maxima 
for d = Api/2 = 182 nm and d — 3Api/2 — 546 nm, which 
indeed agrees very well with the calculated results. Still, 
looking at the behavior of the Green's function is a more 
general way of predicting resonance conditions. 

The results in Fig. [7] calculated with the incident light 
polarized perpendicular to the hole dimer axis show much 
less variation with d. There is a suppression and a red- 
shift of the cross section for d — 40 nm. This is expected 
given the basic behavior illustrated in Fig. Hfb) since in 
this case the field from the induced charges acts to en- 
hance the external field at frequencies below the single- 
hole resonance. But with increasing d the two-hole result 
rather quickly approaches the adjusted one- hole result, 
i.e. the spectrum is only marginally affected by hole-hole 
interactions. This can be anticipated by a look at the 
results for Gyy in Fig. [5] which shows that the long-range 
interaction is rather weak for this configuration. 

Figures [8] and |9] show scattering spectra for chains of 5 
and 8 holes, respectively, illuminated by light polarized 
along the axis of the chain. These results show the same 
trends as those in Fig. [51 but one can still make some 
additional observations, (i) The fact that we have more 
holes means that the collective effects of hole-hole inter- 
actions are stronger since the holes inside the chain now 
have two nearest neighbors. Consequently the peak po- 
sition shifts more now when changing d and the spectra 
rise higher above the (adjusted) one-hole result, (ii) The 
maximum scattering cross section is obtained at some- 
what larger values of d compared with the two-hole case. 
The reason is that not only nearest-neighbor interactions 
matter now. The cross section can be increased by mov- 
ing the next-nearest neighbor hole closer to the second 
maximum of Re[Gxx] at x = 650 nm, see Fig. [5l some- 
thing that is achieved by an increase of d. (iii) We also 
see that the spectral features are sharper here than in 
the two-hole case. This is a rather natural consequence 
of the facts discussed above. An increasing number of 
holes brings an increasing degree of collective behavior 
and constructive interference to the optical response of 
the hole system, which at the same time is more sensi- 
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FIG. 8. (color online) Forward scattering spectra for 5 
nanoholes illuminated by light polarized along the axis of the 
chain of holes. The remaining parameters are the same as in 
Fig. El 



The results presented here agree very well with the 
experimental results found in Ref. [23l see in particular 
Fig. 2 there, (i) As in the experiment nanohole inter- 
actions play an important role when the electric field is 
polarized along the axis of the hole chain, while inter- 
actions only have a minor influence on the spectrum in 
the case of perpendicular polarization, (ii) For parallel 
polarization the experimental scattering spectrum goes 
through the same development as in Figs.[6l[8l and HI For 
the smallest edge-to-edge distances the spectrum is sup- 
pressed and blue-shifted, but as d increases a strong suc- 
cessively red-shifted peak builds up. (iii) The maximum 
scattering cross section in the two-hole case is reached for 
d = 160 nm here, and for d = 150 nm in the experiment, 
see Fig. 2(a) of Ref. [2^. These peak wavelengths differ 
somewhat, ~ 655 nm in the experiment and w 675 nm 
here, part of the reason for this is probably that the holes 
used in the experiment were somewhat smaller, with a 
diameter of 75 nm. (iv) Comparing the experimental re- 
sults for 8 holes with those for 2 holes [Fig 2(a) of the 
experimental paper] we also see much the same trends as 
discussed above. More holes give stronger and sharper 
peaks and bigger wavelength shifts as a function of the 
edge-to-edge distance, just as in the calculation. 



tive to changes in either the photon energy of the incident 
light, the hole-hole separation, or for that matter the di- 
electric environment. Going from 2 holes to 5 makes 
more of a qualitative difference than increasing the num- 
ber from 5 to 8. The reason for this is primarily the 
fact that nearest-neighbor interactions play a dominant 
role; with two holes in the chain both of them just have 
one nearest neighbor, whereas for 5- or 8-hole chains the 
majority of the holes have two nearest neighbors. 
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FIG. 9. (color online) Forward scattering spectra for 8 
nanoholes illuminated by light polarized along the axis of the 
chain of holes. The remaining parameters are the same as in 
Fig.H 



VII. SUMMARY 

In this paper we have presented a derivation of the 
electromagnetic Green's function in systems where the 
background dielectric function varies stepwise along one 
of the coordinate directions, z. The derivation is built on 
a transfer-matrix calculation of the Fourier transform of 
the GF. We have discussed certain symmetry properties 
of the Green's function and also studied its long-range 
properties in real space based on the analytic properties 
of the Fourier transform in the complex plane. 

As an example of an application we have studied the 
long-range properties of the Green's function near a thin 
Au film on a glass substrate. We find there three different 
regimes depending on the lateral distance p between the 
source and field points: (i) A near-field regime where the 
the square of the GF decay as (ii) For 100 nm ^ p ^ 

10 /im the Green's function is dominated by contributions 
from propagating surface plasmons, and |Gp ^ 1/p. (iii) 
Finally, for larger distances, beyond the surface plasmon 
propagation length, the Green's function is dominated 
by contributions from boundary waves (Norton waves) 
grazing the interface. A nearly destructive interference 
between the incident and reflected wave results in the 
intensity cx |Gp decaying as 1/p* in this case. 

We have also applied the Green's function method to a 
calculation of the scattering off two or several nanoholes 
in a thin Au film. We find a strong hole-hole interaction 
mediated by the surface plasmons of the Au film provided 
that the incident electric field is polarized along the axis 
of the hole chain. By increasing the number of holes the 
scattering spectrum gets sharper features and becomes 
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more sensitive to changes in geometry, photon energy, or 
dielectric environment, something that can have apphca- 
tions in for example biochemical sensing. 
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